# The Bombing of Hospitals and Insurgent Responses in Civil War. Evidence from Syria
# Journal of Conflict Resolution

# Regine Schwab (PRIF, Frankfurt/Main), Werner Krause (University of Potsdam), Samer Massoud (Independent researcher)
# 2025-11-26

cat( '\14' )
rm( list = ls( ))

library( dplyr ) # vers. 1.1.4
library( magrittr ) # vers. 2.0.3
library( stringr ) # vers. 1.5.1 
library( purrr ) # vers. 1.0.2
library( PanelMatch ) # vers. 3.1.1
library( ggplot2 ) # vers. 3.5.2

## Variables
## year: Year
## week: Week
## admin2: Admin-2 name
## time_int: Week counter   
## unit_int: Admin2-counter
## hosp_att_bin: Medical facility attack (binary)
## hosp_att_neighb: Medical facility attacks in neighboring unit
## hosp_intent: Intentional medical facility attack (1: intentional)
## reb_vs_reg: Insurgent vs regime military activities
## islglob_vs_reg: Global Islamist vs regime military activities
## reb_vs_reg_fat: Insurgent vs regime military activities (fatalities)
## islglob_vs_reg_fat: Global Islamist vs regime military activities (fatalities)
## reb_vs_reg_fat_log: Insurgent vs regime military activities (fatalities, log( x + 1 ))
## islglob_vs_reg_fat_log: Global Islamist vs regime military activities (fatalities, log( x + 1 ))
## civ_infra: Attacks against civil infrastructure 
## reb_terr_gain: Insurgent territory gains (binary)
## reg_terr_gain: Regime territory gains (binary)
## pop: Population (census)
## unem: Unemployment rate (census)                  
## safe_wat: % of household with safe water (census)
## college: % of college (sensus)
## period: Time period (1: 2017-2018, 0: 2019-2020)
## kurd_fat: Kurdish clashes military activities (fatalities)
## viol_civ: Violence against civilians (binary)             
## ged_reb_vs_reg_fat: Insurgent vs regime military activities (fatalities, GED)
## ged_is_fat: Global Islamist vs regime military activities (fatalities)
## reb_init: Insurgent-initiated military activities
## reb_init_fat: Insurgent-initiated military activities (fatalities)
## reb_init_fat_log: Insurgent-initiated military activities (fatalities, log( x + 1 ))
## reb_init_wmines_fat: Insurgent-initiated military activities with explosives (fatalities)
## reb_init_wterr_fat: Insurgent-initiated military activities with territory gains (fatalities)
## reg_init: Regime-initiated military activities
## reg_init_fat: Regime-initiated military activities (fatalities)
## reg_init_fat_log: Regime-initiated military activities (fatalities, log( x + 1 ))
## reg_init_wterr_fat: Regime-initiated military activities with territory gains (fatalities)
## ratio_fat: Ratio of insurgent- and regime-initiated activities (fatalities)
## unid_init: Military activities with unclear initiator
## unid_init_fat: Military activities with unclear initiator (fatalities)
## unid_init_fat_log: Military activities with unclear initiator (fatalities, log( x + 1 ))
## unid_init_woterr_fat: Military activities with unclear initiator w/o territory gains (fatalities)
## geo_prec_mean: Average precision of location of military activities
## geo_prec_sh3: Share of military activities with imprecise locations (ACLED-code: 3)
## geo_prec_sh23: Share of military activities with imprecise locations (ACLED-code: 2 & 3)
## time_prec_mean: Average precision of time of military activities
## time_prec_sh3: Share of military activities with imprecise time (ACLED-code: 3)
## time_prec_sh23: Share of military activities with imprecise time (ACLED-code: 2 & 3)

## Note: Due to the sensitivity of the data, Figures 2[B], A8, and A14 are not included. 
## Please contact authors for additional information.

folder.path <- '' # <- set your path
d <- readRDS( paste0( folder.path , '_bhircw.rds' ))

outc <- 'reb_vs_reg_fat'
treat <- 'hosp_att_bin'
covs <- c( 'civ_infra' , 'islglob_vs_reg_fat' , 'reb_terr_gain' , 'reg_terr_gain' , 'hosp_att_neighb' , 'pop' , 'unem' , 'safe_wat' , 'college' )

c.formula <- paste0( '~ I( lag( ' , outc , ', 1:4)) + ' ,  paste0( 'I( lag( ' , covs , ', 1:4))' , collapse = ' + '))


PM_wrap <- function( data , lag = 4 , lead = 0 : 4 , refinement.method = 'ps.weight' , covs.formula , qoi = 'att' 
                     , forbid.treatment.reversal = F , match.missing = F , moderator = NULL ){
  out <- list( ) 
  out$PM <- PanelMatch( panel.data = data , lag = lag , lead = lead , refinement.method = refinement.method
                        , covs.formula = covs.formula , qoi = qoi , forbid.treatment.reversal = forbid.treatment.reversal 
                        , match.missing = match.missing )
  set.seed( 78208 )
  out$PE$CI95 <- PanelEstimate( out$PM , panel.data = data 
                                , moderator = moderator 
                                , confidence.level = .95 )
  set.seed( 78208 )
  out$PE$CI90 <- PanelEstimate( out$PM , panel.data = data 
                                , moderator = moderator 
                                , confidence.level = .9 )
  return( out )
}

# Table A2 ====

d$weeks %>%
  select( any_of( outc ) , any_of( covs ) , 'reb_init_fat' , 'reg_init_fat' , 'unid_init_fat'  ) %>%
  na.omit( ) %>%
  psych::describe( )

# Figure 2 ====

d$weeks %>%
  ggplot( aes( x = time_int , y = admin2 , fill = factor( hosp_att_bin ) , color = factor( hosp_att_bin ) )) +
  geom_tile( ) +
  scale_y_discrete(limits=rev) +
  scale_fill_manual( values = c( '#a0a0a001' , 'grey50' )) +
  scale_color_manual( values = c( '#a0a0a001' , 'black' )) 

# Figure 3 ===== 

o <- PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
              , covs.formula = as.formula( c.formula ))

summary( o$PE$CI95$matched.sets )$overview %>% 
  ggplot( aes( x = matched.set.size )) +
  geom_histogram( binwidth = 1 , fill = 'grey70' , color = 'black') +  
  scale_x_continuous( breaks = seq( 0 , 17 , 1 ) , name = '\nNumber of Matched Control Units' ) + 
  scale_y_continuous( name = 'Frequency\n' , breaks = seq( 0 , 12 , 1 )) 

# Figure 4 =====

cov_balance <- get_covariate_balance( o$PM
                                      , panel.data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                                      , covariates = c( outc , covs )) %>% 
  as.data.frame( ) %>% .[ c(1 : 4) , ] %>%
  mutate( type = 'ps.weight' 
          , time = row.names( . )
          , time = stringr::str_remove( time , 't_' )
          , time = ifelse( time > 0 , as.numeric( time ) * -1 , as.numeric( time ))) %>%
  tidyr::gather( var , val , all_of( paste0( 'o.PM.att.' , outc )) , all_of( paste0( 'o.PM.att.' , covs )))

cov_balance %>%
  ggplot( aes( time , val )) +
  geom_line( size = 1 ) +
  facet_wrap( ~var , nrow = 2 ) +
  scale_y_continuous( name = 'Standardized mean difference\n' , limits = c( -1.5 , 1.5 ) ) +
  scale_x_continuous( name = '\nWeeks prior to treatment' ) +
  geom_hline( yintercept = .25 ) +
  geom_hline( yintercept = -.25 )

# Figure 5 ===== 

summary( o$PE$CI95 ) %>% as.data.frame( ) %>%
  mutate( time = row.names( . )) %>%
  ggplot( aes( x = time , y = estimate )) +
  geom_point( ) +
  geom_errorbar( aes( ymin = `2.5%` , ymax = `97.5%` ) , width = 0 )

# Figure A1 ====

o <- map_df( c( 'none' , 'mahalanobis' , 'ps.match' , 'CBPS.match' , 'ps.weight' , 'CBPS.weight' )
             , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                          , covs.formula = as.formula( c.formula ) , refinement.method = .x )$PM %>% 
               get_covariate_balance( . , panel.data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                                      , covariates = c( outc , covs )) %>% 
               as.data.frame( ) %>% .[ c(1 : 4) , ] %>%
               mutate( type = .x , time = row.names( . ) , time = stringr::str_remove( time , 't_' )
                       , time = ifelse( time > 0 , as.numeric( time ) * -1 , as.numeric( time ))) %>%
               tidyr::gather( var , val , all_of( paste0( '..att.' , outc )) , all_of( paste0( '..att.' , covs ))))

o %>%
  filter( time != 0 ) %>%
  ggplot( aes( time , val , color = var )) +
  geom_line( ) +
  facet_wrap( ~type ) +
  scale_y_continuous( name = 'Standardized mean difference\n' , limits = c( -1 , 1 ) ) +
  scale_x_continuous( name = '\nWeeks prior to treatment' ) +
  geom_hline( yintercept = .25 ) +
  geom_hline( yintercept = -.25 ) +
  theme( legend.position = 'bottom' )

# Figure A2 ====

o <- map_df( c( 'none' , 'mahalanobis' , 'ps.match' , 'CBPS.match' , 'ps.weight' , 'CBPS.weight' )
             , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                          , covs.formula = as.formula( c.formula ) , refinement.method = .x )$PE$CI95 %>% 
               summary( . ) %>%
               as.data.frame( ) %>%
               mutate( type = .x
                       , time = row.names( . )))

o %>% 
  ggplot( aes( x = time , y = estimate  , group = type , color = type )) +
  geom_point( position = position_dodge( width = 0.5 ) , size = 1.75 ) +
  geom_errorbar( aes( ymin = `2.5%` , ymax = `97.5%` ) , position = position_dodge( width = 0.5 ) , width = 0 ) +
  geom_hline( yintercept = 0 )

# Figure A3 ====

o <- map_df( seq( 1 , 12 , 1 )
             , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                          , covs.formula = as.formula( c.formula ) , lag = .x )$PE$CI95 %>% 
               summary( . ) %>%
               as.data.frame( ) %>%
               mutate( type = .x
                       , time = row.names( . )))

o %>%
  ggplot( aes( x = time , y = estimate  , group = type , color = type )) +
  geom_point( position = position_dodge( width = 0.5 ) , size = 1.75 ) +
  geom_errorbar( aes( ymin = `2.5%` , ymax = `97.5%` ) , position = position_dodge( width = 0.5 ) , width = 0 ) +
  geom_hline( yintercept = 0 )

# Figure A4 ====

o <- map_df( unique( d$weeks$admin2 )
             , ~ PM_wrap( data = d$weeks %>% filter( admin2 != .x ) %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                          , covs.formula = as.formula( c.formula ))$PE$CI95 %>% 
               summary( . ) %>%
               as.data.frame( ) %>%
               mutate( type = .x
                       , time = row.names( . )))

o %>% 
  ggplot( aes( x = time , y = estimate  , group = type , color = type )) +
  geom_point( position = position_dodge( width = 0.5 ) , size = 1.75 ) +
  geom_errorbar( aes( ymin = `2.5%` , ymax = `97.5%` ) , position = position_dodge( width = 0.5 ) , width = 0 ) +
  geom_hline( yintercept = 0 )

# Figure A5 ====

o <- map_df( c( 0 , 1 )
             , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                          , covs.formula = as.formula( c.formula )
                          , moderator = "period" )$PE$CI95[[ as.character( .x )]] %>%
               summary( ) %>%
               as.data.frame( ) %>%
               mutate( time = row.names( . ) , type = .x )) %>%
  mutate( type = as.factor( type ))

o %>% 
  ggplot( aes( x = time , y = estimate  , group = type , color = type )) +
  geom_point( position = position_dodge( width = 0.5 ) , size = 1.75 ) +
  geom_errorbar( aes( ymin = `2.5%` , ymax = `97.5%` ) , position = position_dodge( width = 0.5 ) , width = 0 ) +
  geom_hline( yintercept = 0 )

# Figure A6 ====

o <- map_df( c( 0 , 1 )
             , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                          , covs.formula = as.formula( c.formula )
                          , moderator = "hosp_intent" )$PE$CI95[[ as.character( .x )]] %>%
               summary( ) %>%
               as.data.frame( ) %>%
               mutate( time = row.names( . ) , type = .x )) %>%
  mutate( type = as.factor( type ))

o %>% 
  ggplot( aes( x = time , y = estimate  , group = type , color = type )) +
  geom_point( position = position_dodge( width = 0.5 ) , size = 1.75 ) +
  geom_errorbar( aes( ymin = `2.5%` , ymax = `97.5%` ) , position = position_dodge( width = 0.5 ) , width = 0 ) +
  geom_hline( yintercept = 0 )

# Table A3 ====

## Model 1
PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , 'reb_vs_reg_fat_log' ) 
         , covs.formula = as.formula( str_replace_all( c.formula , '_fat' , '_fat_log' )))$PE$CI95 %>% summary( )

## Model 2
PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , 'reb_vs_reg' ) 
         , covs.formula = as.formula( str_remove_all( c.formula , '_fat' )))$PE$CI95 %>% summary( )

## Model 3
PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc ) 
         , covs.formula = as.formula( str_remove( c.formula , ' I\\( lag\\( reb_vs_reg_fat, 1:4\\)\\) \\+' )))$PE$CI95 %>% summary( )

# Figure A7 ====

o <- map_df( c( 'geo_prec_mean' , 'geo_prec_sh3' , 'geo_prec_sh23'
                , 'time_prec_mean' , 'time_prec_sh3' , 'time_prec_sh23' )
             , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , .x )
                          , covs.formula = as.formula( paste0( c.formula , '+ I( lag( ' , .x , ', 1:4 ))'
                          ))
                          , match.missing = TRUE )$PE$CI95 %>%
               summary( ) %>%
               as.data.frame( ) %>%
               mutate( time = row.names( . ) , type = .x ))

o %>% 
  ggplot( aes( x = time , y = estimate )) +
  geom_point( position = position_dodge( width = 0.5 ) , size = 1.75 ) +
  geom_errorbar( aes( ymin = `2.5%` , ymax = `97.5%` ) , position = position_dodge( width = 0.5 ) , width = 0 ) +
  geom_hline( yintercept = 0 ) + 
  facet_wrap( ~ type  )

# Table A4 ====

## Model 1
PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , 'islglob_vs_reg_fat' ) 
         , covs.formula = as.formula( c.formula ))$PE$CI95 %>% summary( )

## Model 2
PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , 'kurd_fat' ) 
         , covs.formula = as.formula( paste0( c.formula , ' + I( lag( kurd_fat, 1:4))' )))$PE$CI95 %>% summary( )

## Model 3
PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , 'viol_civ' ) 
         , covs.formula = as.formula( paste0( c.formula , ' + I( lag( viol_civ, 1:4))' )))$PE$CI95 %>% summary( )

# Table A5 ====

## Model 1
c.temp <- c( 'ged_reb_vs_reg_fat' , 'ged_is_fat' , 'hosp_att_neighb' , 'pop' , 'unem' , 'safe_wat' , 'college' )
c.temp <- paste0( '~ ' ,  paste0( 'I( lag( ' , c.temp , ', 1:4 ))' , collapse = ' + ' ))
PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , 'ged_reb_vs_reg_fat' ) 
         , covs.formula = as.formula( c.temp ))$PE$CI95 %>% summary( )

## Model 2
PM_wrap( data = d$weeks %>% filter( max( hosp_att_bin ) == 1 , .by = admin2 ) %>% PanelData( . , 'unit_int' , 'time_int' , treat , 'ged_reb_vs_reg_fat' ) 
         , covs.formula = as.formula( c.temp ))$PE$CI95 %>% summary( )

## Model 3
PM_wrap( data = d$days %>% PanelData( . , 'unit_int' , 'time_int' , treat , 'ged_reb_vs_reg_fat' ) 
         , covs.formula = ~ I( lag( ged_reb_vs_reg_fat , 1:14 )) + 
           I( lag( ged_is_fat , 1:14 )) +
           I( lag( hosp_att_neighb , 1:14 )) +
           I( lag( pop, 1:14 )) + 
           I( lag( unem, 1:14 )) + 
           I( lag( safe_wat, 1:14 )) + 
           I( lag( college, 1:14 ))
         , lag = 14
         , lead = 0 : 14 )$PE$CI95 %>% summary( )

# Table A6 ====

## Model 1
PM_wrap( data = d$weeks %>% filter( max( hosp_att_bin ) == 1 , .by = admin2 ) %>% 
           PanelData( . , 'unit_int' , 'time_int' , treat , outc ) 
         , covs.formula = as.formula( c.formula ))$PE$CI95 %>% summary( )

## Model 2
PM_wrap( data = d$days %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc ) 
         , covs.formula = ~ I( lag( reb_vs_reg_fat , 1:14 )) + 
           I( lag( civ_infra , 1:14 )) + 
           I( lag( reb_terr_gain , 1:14 )) + 
           I( lag( reg_terr_gain , 1:14 )) + 
           I( lag( islglob_vs_reg_fat , 1:14 )) + 
           I( lag( hosp_att_neighb , 1:14 )) +
           I( lag( pop, 1:14 )) + 
           I( lag( unem, 1:14 )) + 
           I( lag( safe_wat, 1:14 )) + 
           I( lag( college, 1:14 ))
         , lag = 14
         , lead = 0 : 14 )$PE$CI95 %>% summary( )

## Model 3-6
c.temp <- c( 'civ_infra' , 'islglob_vs_reg_fat' , 'reb_terr_gain' , 'reg_terr_gain' , 'hosp_att_neighb' )
c.formula.temp <- paste0( '~ I( lag( ' , outc , ', 1:4)) + ' ,  paste0( 'I( lag( ' , c.temp , ', 1:4))' , collapse = ' + '))

## Model 3
PM_wrap( data = d$grids[[ 1 ]] %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc ) , covs.formula = as.formula( c.formula.temp ))$PE$CI95 %>% summary( )

## Model 4
PM_wrap( data = d$grids[[ 2 ]] %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc ) , covs.formula = as.formula( c.formula.temp ))$PE$CI95 %>% summary( )

## Model 5
PM_wrap( data = d$grids[[ 3 ]] %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc ) , covs.formula = as.formula( c.formula.temp ))$PE$CI95 %>% summary( )

## Model 6
PM_wrap( data = d$grids[[ 4 ]] %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc ) , covs.formula = as.formula( c.formula.temp ))$PE$CI95 %>% summary( )

# Table 1 =====

covs <- c( 'reb_init_fat' , 'reg_init_fat' , 'unid_init_fat' , 'civ_infra' , 'islglob_vs_reg_fat' , 'reb_terr_gain' , 'reg_terr_gain' , 'hosp_att_neighb'
           , 'pop' , 'unem' , 'safe_wat' , 'college' )
c.formula <- paste0( '~ ' ,  paste0( 'I( lag( ' , covs , ', 1:4 ))' , collapse = ' + '))
map( c( 'reb_init_fat' , 'reg_init_fat' , 'unid_init_fat' )
     , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , .x ) , covs.formula = as.formula( c.formula ) )$PE$CI95 %>% 
       summary( . )
)

# Figure A9 =====

outc <- 'reb_init_fat'
o <- map_df( c( 'none' , 'mahalanobis' , 'ps.match' , 'CBPS.match' , 'ps.weight' , 'CBPS.weight' )
             , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                          , covs.formula = as.formula( c.formula ) , refinement.method = .x )$PM %>% 
               get_covariate_balance( . , panel.data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                                      , covariates = c( outc , covs )) %>% 
               as.data.frame( ) %>% .[ c(1 : 4) , ] %>%
               mutate( type = .x , time = row.names( . ) , time = stringr::str_remove( time , 't_' )
                       , time = ifelse( time > 0 , as.numeric( time ) * -1 , as.numeric( time ))) %>%
               tidyr::gather( var , val , all_of( paste0( '..att.' , outc )) , all_of( paste0( '..att.' , covs ))))

o %>%
  filter( time != 0 ) %>%
  ggplot( aes( time , val , color = var )) +
  geom_line( ) +
  facet_wrap( ~type ) +
  scale_y_continuous( name = 'Standardized mean difference\n' , limits = c( -1 , 1 ) ) +
  scale_x_continuous( name = '\nWeeks prior to treatment' ) +
  geom_hline( yintercept = .25 ) +
  geom_hline( yintercept = -.25 ) +
  theme( legend.position = 'bottom' )

# Figure A10 =====

o <- map_df( c( 'none' , 'mahalanobis' , 'ps.match' , 'CBPS.match' , 'ps.weight' , 'CBPS.weight' )
             , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                          , covs.formula = as.formula( c.formula ) , refinement.method = .x )$PE$CI95 %>% 
               summary( . ) %>%
               as.data.frame( ) %>%
               mutate( type = .x , time = row.names( . )))

o %>% 
  ggplot( aes( x = time , y = estimate  , group = type , color = type )) +
  geom_point( position = position_dodge( width = 0.5 ) , size = 1.75 ) +
  geom_errorbar( aes( ymin = `2.5%` , ymax = `97.5%` ) , position = position_dodge( width = 0.5 ) , width = 0 ) +
  geom_hline( yintercept = 0 )

# Figure A11 ====

o <- map_df( seq( 1 , 12 , 1 )
             , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                          , covs.formula = as.formula( c.formula ) , lag = .x )$PE$CI95 %>% 
               summary( . ) %>%
               as.data.frame( ) %>%
               mutate( type = .x , time = row.names( . )))

o %>% 
  ggplot( aes( x = time , y = estimate  , group = type , color = type )) +
  geom_point( position = position_dodge( width = 0.5 ) , size = 1.75 ) +
  geom_errorbar( aes( ymin = `2.5%` , ymax = `97.5%` ) , position = position_dodge( width = 0.5 ) , width = 0 ) +
  geom_hline( yintercept = 0 )

# Figure A12 ====

o <- map_df( unique( d$weeks$admin2 )
             , ~ PM_wrap( data = d$weeks %>% filter( admin2 != .x ) %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                          , covs.formula = as.formula( c.formula ))$PE$CI95 %>% 
               summary( . ) %>%
               as.data.frame( ) %>%
               mutate( type = .x , time = row.names( . )))

o %>% 
  ggplot( aes( x = time , y = estimate  , group = type , color = type )) +
  geom_point( position = position_dodge( width = 0.5 ) , size = 1.75 ) +
  geom_errorbar( aes( ymin = `2.5%` , ymax = `97.5%` ) , position = position_dodge( width = 0.5 ) , width = 0 ) +
  geom_hline( yintercept = 0 )

# Figure A13 ====

o <- map_df( c( 0 , 1 )
             , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                          , covs.formula = as.formula( c.formula )
                          , moderator = "period" )$PE$CI95[[ as.character( .x )]] %>%
               summary( ) %>% as.data.frame( ) %>%
               mutate( time = row.names( . ) , type = .x )) %>%
  mutate( type = as.factor( type ))

o %>% 
  ggplot( aes( x = time , y = estimate  , group = type , color = type )) +
  geom_point( position = position_dodge( width = 0.5 ) , size = 1.75 ) +
  geom_errorbar( aes( ymin = `2.5%` , ymax = `97.5%` ) , position = position_dodge( width = 0.5 ) , width = 0 ) +
  geom_hline( yintercept = 0 )

# Figure A14 ==== 

o <- map_df( c( 0 , 1 )
             , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc )
                          , covs.formula = as.formula( c.formula )
                          , moderator = "hosp_intent" )$PE$CI95[[ as.character( .x )]] %>%
               summary( ) %>% as.data.frame( ) %>%
               mutate( time = row.names( . ) , type = .x )) %>%
  mutate( type = as.factor( type ))

o %>% 
  ggplot( aes( x = time , y = estimate  , group = type , color = type )) +
  geom_point( position = position_dodge( width = 0.5 ) , size = 1.75 ) +
  geom_errorbar( aes( ymin = `2.5%` , ymax = `97.5%` ) , position = position_dodge( width = 0.5 ) , width = 0 ) +
  geom_hline( yintercept = 0 )

# Table A7 ====

## Model 1
PM_wrap( data = d$weeks %>% filter( max( hosp_att_bin ) == 1 , .by = unit_int ) %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc ) 
         , covs.formula = as.formula( c.formula ) )$PE$CI95 %>% summary( )

## Model 2
PM_wrap( data = d$days %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc ) 
         , covs.formula = ~ I( lag( reb_vs_reg_fat , 1:14 )) +
           I( lag( islglob_vs_reg_fat , 1:14 )) +
           I( lag( hosp_att_neighb , 1:14 )) +
           I( lag( pop, 1:14 )) + 
           I( lag( unem, 1:14 )) + 
           I( lag( safe_wat, 1:14 )) + 
           I( lag( college, 1:14 ))
         , lag = 14
         , lead = 0 : 14 )$PE$CI95 %>% summary( )

## Model 3-6
covs.temp <- c( 'reb_init_fat' , 'reg_init_fat' , 'unid_init_fat' , 'civ_infra' , 'islglob_vs_reg_fat' , 'reb_terr_gain' , 'reg_terr_gain' , 'hosp_att_neighb' )
c.formula.temp <- paste0( '~ ' ,  paste0( 'I( lag( ' , covs.temp , ', 1:4 ))' , collapse = ' + '))

## Model 3
PM_wrap( data = d$grids[[ 1 ]] %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc ) , covs.formula = as.formula( c.formula.temp ))$PE$CI95 %>% summary( )

## Model 4
PM_wrap( data = d$grids[[ 2 ]] %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc ) , covs.formula = as.formula( c.formula.temp ))$PE$CI95 %>% summary( )

## Model 5
PM_wrap( data = d$grids[[ 3 ]] %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc ) , covs.formula = as.formula( c.formula.temp ))$PE$CI95 %>% summary( )

## Model 6
PM_wrap( data = d$grids[[ 4 ]] %>% PanelData( . , 'unit_int' , 'time_int' , treat , outc ) , covs.formula = as.formula( c.formula.temp ))$PE$CI95 %>% summary( )

# Table A8 ====

PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , 'reb_init_wmines_fat' ) , covs.formula = as.formula( c.formula ) )$PE$CI95 %>% summary( )

# Table A9 ====

map( c( 'reb_init_wterr_fat' , 'reg_init_wterr_fat' , 'unid_init_woterr_fat' )
     , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , .x ) , covs.formula = as.formula( c.formula ) )$PE$CI95 %>% 
       summary( . ))

# Table A10 ====

map( c( 'reb_init_fat_log' , 'reg_init_fat_log' , 'unid_init_fat_log' )
     , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , .x ) 
                  , covs.formula = as.formula( str_replace_all( c.formula , '_fat' , '_fat_log' )))$PE$CI95 %>% 
       summary( . ))

# Table A11 ====

map( c( 'reb_init' , 'reg_init' , 'unid_init' )
     , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , .x ) 
                  , covs.formula = as.formula( str_remove_all( c.formula , '_fat' )))$PE$CI95 %>% 
       summary( . ))

# Table A12 ====

map( c( 'reb_init_fat' , 'reg_init_fat' , 'unid_init_fat' )
     , ~ PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , .x ) 
                  , covs.formula = as.formula( str_remove_all( c.formula , paste0( ' I\\( lag\\( ' , .x , ', 1:4\\)\\) \\+' ))))$PE$CI95 %>% 
       summary( . ))

# Table A13 ====

covs.temp <- c( 'unid_init_fat' , 'civ_infra' , 'islglob_vs_reg_fat' , 'reb_terr_gain' , 'reg_terr_gain' , 'hosp_att_neighb' 
                , 'pop' , 'unem' , 'safe_wat' , 'college' )
c.formula.temp <- paste0( '~ ' ,  paste0( 'I( lag( ' , covs.temp , ', 1:4 ))' , collapse = ' + '))
PM_wrap( data = d$weeks %>% PanelData( . , 'unit_int' , 'time_int' , treat , 'ratio_fat' ) , covs.formula = as.formula( c.formula.temp ) )$PE$CI95 %>% summary( )
